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Abstract. We report on a comparison of high-resolution numerical simulations 
of Lagrangian particles advected by incompressible turbulent hydro- and magne- 
tohydrodynamic (MHD) flows. Numerical simulations were performed with up to 
1024^ collocation points and 10 million particles in the Navier-Stokes case and 512'^ 
collocation points and 1 million particles in the MHD case. In the hydrodynamics 
case our findings compare with recent experiments from Mordant et al. [1] and 
Xu et al. [2]. They differ from the simulations of Biferale et al. [3] due to differ- 
ences of the ranges choosen for evaluating the structure functions. In Navier-Stokes 
turbulence intermittency is stronger than predicted by a multifractal approach of 
[3] whereas in MHD turbulence the predictions from the multifractal approach are 
more intermittent than observed in our simulations. In addition, our simulations 
reveal that Lagrangian Navier-Stokes turbulence is more intermittent than MHD 
turbulence, whereas the situation is reversed in the Eulerian case. Those findings 
can not consistently be described by the multifractal modeling. The crucial point 
is that the geometry of the dissipative structures have different implications for 
Lagrangian and Eulerian intermittency. Application of the multifractal approach 
for the modeling of the acceleration PDFs works well for the Navier-Stokes case 
but in the MHD case just the tails are well described. 



1. Introduction 

Lagrangian statistics of turbulent flows has undergone a rapid development in the 
last 6 years due to enormous progress in experimental techniques measuring parti- 
cle trajectories. Particle tracking velocimetry has been used for moderate Reynolds 
numbers by [4]. However, the techniques developed in Cornell [5, 6, 7] and Lyon [8, 1] 
allowed the measurements of probability density functions (PDFs) of velocity in- 
crements which triggered a renewed interest in the theoretical understanding of 
Lagrangian statistics. A very promising approach based on a Markovian closure 
was recently introduced in [9] . This approach is not readily applicable to magneto- 
hydrodynamic (MHD) turbulence, since the distribution function depends not only 
on velocity, space and initial condition but in addition on the Jacobian. Although 
work in this direction is in progress, we compare our simulations to a phenomeno- 
logical model of Lagrangian statistics in Navier-Stokes turbulence introduced by 
Biferale et al. [3]. Our findings show increased intermittency in Navier-Stokes flows. 
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such that the structure functions agree with recent experimental data from the two 
experimental groups [1, 2]. On the other hand, our MHD simulations are less inter- 
mittent than the predictions from a multifractal model. 



2. Numerical Methods 

The Lagrangian particle trajectories were obtained by two slightly different parallel 
spectral codes (Garching and Bochum) based on the spectral code used in [10]. 
The velocity and magnetic field was evaluated at the particle positions using either 
trilinear or tricubic interpolation. Contrary to the simulations of [3] , we found that 
tricubic interpolation captures especially trajectories with high acceleration more 
precise than linear interpolation, a conclusion which was drawn already 20 years 
ago (see [11] and the discussion therein). However, the effect of the slightly different 
trajectories has only a minor effect on the tails of the acceleration PDFs with the 
tendency that the PDF calculated with tricubic interpolation is slightly more in- 
termittent than the corresponding one calculated with trilinear interpolation. The 
simulations presented here use a tricubic interpolation for obtaining the velocities 
at the particle positions. 



3. Navier-Stokes turbulence 

We performed a set of simiilations for the Navier-Stokes-equations with two reso- 
lutions, 512^ and 1024^ collocation points and 1 million and 10 million particles, 
respectively, in order to obtain reliable statistical results within a few large eddy- 
turnover times. 

Parameters of all simulations (Navier-Stokes and MHD) arc summarized in Ta- 
ble 1. Here we used the same conventions as described in [12]. To get to a stationary 
turbulent state, we started with randomly distributed Fourier-modes without driv- 
ing. After the turbulence has been fully developed, the low mode number modes 
were kept constant. Particles with initially homogeneous random positions were 
injected when a stationary state was reached. For the simulation with 1024^ collo- 
cation points, we started from the stationary turbulent state with 512'^ collocation 
points. The relaxation to a new stationary state took about one largc-cddy turnover 
time. After this period, the particles were injected. 

The choice of the parameters for Runl, Run4 and Run5 and therefore the result- 
ing Taylor microscale Reynolds number was motivated by the standard procedure 
to choose the dissipation length Id smaller than the grid spacing [13]. However, 
recent investigations of how dissipative structures like shocks, tubes and sheets en- 
ter the dissipation range [14, 15, 16] suggest a more conservative choice which was 
realized in Run2 and Run3. 

The exponents of the longitudinal Eulerian structure functions Sp =< |u(x-|-l) — 
u(x)|P >, angular brackets denoting spatial averaging, can be described by the She- 
Lcveque formula [17, 18] for the simulations Runl - Rim3. The structure functions 
using extended self similarity (ESS) [19] together with straight lines illustrating the 
fitted slopes are shown in Figure 1 for data obtained from Run3. The values of 
the corresponding exponents are summarized in Table 2. These Eulerian statistics 
serve just as a test to check the numerics. The determination of the Lagrangian 
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Table 1. Parameters of the numerical simulations. Rx: Taylor microscale Reynolds number 
y^^TSu^L/u, uo = ^y2/3Ek, Ek- kinetic energy, Em- magnetic energy, E ^ Ek + Em, ^k- 
kinetic energy dissipation rate, e™: magnetic energy dissipation rate, e = ek + em, v: 
viscosity, 77: resistivity, la'- dissipation lengthscale (i'^/efc)i/i, r^: Kolmogorov time scale 
(^/efe)i/2, L = (2/3E)^^^ /e: integral scale, Tl = L/uq- large-eddy turnover time, T: total 
integration time, N^: number of collocation points, Np-. number of particles, Navier-Stokes 
simulations: Runl-Run3, MHD simulations: Run4, Run5 
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Figure 1. Eulerian structure functions of order p=l,2,4,5,6 (top to bottom) 
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Table 2. Eulerian structure functions obtained using ESS with 1^3 = 1. 
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Figure 2. 2nd order structure function normalized to er for Runl (solid line) and Run4 

(dashed line) 
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Figure 3. Lagrangian structure functions from Run3 of order p=3-6 (top to bottom) 

structure functions Sp =< \u{t + T) — u(t)|'° >, angular brackets denoting temporal 
averaging, turned out to be much more difficult. Figure 2 shows a typical plot 
of the second order Lagrangian structure function normalized to er. It is clear 
that no scaling range is present as already observed in the experiments [1, 2] and 
simulation [3, 20]. Therefore, in order to obtain scaling exponents one has to rely 
on the assumption of extended self similarity. Figure 3 shows an evaluation of the 
Lagrangian structure functions assuming ESS. In Table 3 we present the relative 
exponents Cp/C'i for our simulations. In addition, this table contains also exponents 
collected from the experiments [1, 2] and other numerical simulations [3, 1]. 

First, we observe that the simulations Runl - Run3 all give exponents which 
agree within the error bars. Thus, a possible dependence of the exponents on the 
Reynolds number or the choice of a stricter criterion for the numerical resolution 
could not be detected given the relative large error bars. In addition, the exponents 
fit quite well to the present experiments [1, 2] but are clearly different from the 
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Runl 
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Experiment [2] 



0.58±0.006 
0.57±0.007 
0.57±0.005 
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1.76±0.11 
2.26 
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1.9±0.2 
1.69±0.2 
1.66±0.53 



Table 3. Relative ESS-exponents calculated with respect to the structure function of 
order 2. MF-NS denotes the Multi- Fractal approach for Navier-Stokes. Reynolds numbers 
for the different simulations and experiments: Runl: Rx = 190, Run2: R\ = 122, Run3: 
Rx = 178, Simulation [3]: Rx = 284, Experiment [1]: Rx = 1000, Simulation [1]: Rx = 140, 
Experiment [2]: Rx = 815 



exponents obtained by [3] . The reason for this is not the different interpolation for 
particle velocities (trilinear in [3], tricubic here). We repeated a simulation with 
the parameters of Run2 but using trilinear interpolation and got the same scaling 
as with tricubic interpolation. We explain this discrepancy similar as in [2] by ob- 
serving that the evaluation of the structure functions was performed systematically 
using larger values of r compared to us. We have have chosen the inertial range 
from Figure 2 by the requirement that the function stays above ninety percent of its 
maximum value. This leads to an inertial range of 2 < r < 7. We are able to repro- 
duce the exponents of [3] if we choose a minimal value of r of about eight. Larger 
values of r result in a more Gaussian behavior and explains why the exponents of 
[3] are less intermittent. 

We also applied the multifractal model to the acceleration statistics obtained from 
the simulations. We shortly review the approach of [3] . One starts with a suitable 
description for the Eulerian structure functions, e.g. the She-Leveque model [17] and 
performs a Legendre transformation to obtain the singularity spectrum. In order 
to translate Eulerian to Lagrangian increments, one assumes a Kolmogorov like 
relation SrV ^ Siu where temporal and spatial increments are related by r; ~ l/Siu. 
The resulting expression for the Lagrangian structure functions, 

hp + 3-D(h) 

. (3-1) 

is evaluated by a saddle point integration. To obtain the acceleration PDF, first 
the acceleration is defined as a = where = r^(/i, wq) is the Kolmogorov time 

scale which is itself a multifractal quantity depending on the large scale velocity 
field Uq. Assuming a Gaussian statistics of Uq and integrating over the possible 
scaling factors h results in an explicit expression for the acceleration PDF (see [3] 
for details) 

P{a) ^ dh a^(h-6+D(h))/3)j^y{h) 

exp(^-ia2(i+/.)/3^.(^)^ 



with a = a/aa, (Ja = {a^Y'^ , y{h) = x{h - 5 + D{h))/6 + 2{2D{h) + 2h - 7)/3 
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and z{h) — x(l + ft-)/3 + A{2h — l)/3. A comparison of the multifractal prediction to 
the numericaUy obtained acceleration PDFs is shown in Figure 4. Ahhough there 
is an excehent agreement between prediction and simulation, one has to keep in 
mind that the multifractal prediction contains three parameters. The first param- 




10 20 30 40 50 

a/Og 



Figure 4. Acceleration PDFs for runs Run2(+) and Run3(x) including the PDFs using 
the multifractal approach for the Navier-Stokes simulations. The inset shows the PDF 
multiplied by (a/a)"'. 

eter is hidden in the relation ^a^) cx i?^ when normalizing the width of the PDF. 
However, this parameter can be determined from one simulation and should then 
be kept fixed for other Reynolds numbers. In the She-Leveque model, the value of 
/imin is given by /imin — 1/9. If one uses this value, it is not possible to get good 
agreement with the measured shape of the PDF. Therefore, as in [3], we use /imin 
as a free parameter. The last is a free amplitude in the normalization. In order to 
get such an excellent agreement in Figure 4, we had to choose /imin = 0.175 for 
Run2 and /imin = 0.16 for Rmi3. The dependence on /imax is negligible. 



4. MHD turbulence 

The parameters for the MHD simulations are summarized in Table 1 (Run4, Run5). 
Both runs were performed with negligible magnetic and cross helicity. lESS plots 
of the Lagrangian velocity structure functions are shown in Figure 5. They show 
a similar curved shape although no trapping in vortex tubes appears in MHD 
turbulence (see also Biferale et al. [21]). 

The exponents for the Lagrangian velocity structure functions are given in Ta- 
ble 4. Also shown is the prediction by a multifractal model, which was obtained 
using the same steps as described above, but starting with a She-Leveque like for- 
mula suitable for incompressible MHD-turbulence [22, 10], 

cr°b)-| + i-Q)'^' . (4.1) 

Although this formula is strictly valid only for the structure functions of the Elsasser 
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Figure 5. Lagrangian structure functions from Run4 of order 1,3,4,5,6 (top to bottom) 

variables — u ± B, we assume a cascade in the kinetic energy so that this 
formula can also be applied to the structure functions of velocity. The resulting 
multifractal model shows now an increased degree of intermittency compared to 
the numerical simulations (Run4, Run5). On the first sight this is astonishing since 
this is just the opposite behavior as in the Navier-Stokes case. To summarize, we 
have the following situation that in the Eulerian description, MHD turbulence is 
more intermittent than Navier-Stokes turbulence whereas the situation is reversed 
in the Lagrangian picture. This finding is also not compatible with the multifractal 
ansatz. The multifractal ansatz possesses a certain monotonicity property. This 
means that if for two different sets of structure function exponents, one is more 
intermittent than the other in the Eulerian picture, than this one is also more 
intermittent in the Lagrangian turbulence. To see this, it is sufficient to look at 
high values p of the order of the structure functions. One observes that the value 
of h* where the infimum of 

hp + 2,- D{h) 

is assumed goes to /imin for high values of p. Thus the asymptotic behavior reads 

Cp = /iminP + 3 - D{hrain) , p > 1 . 

For the saddle point evaluation of the Lagrangian structure functions (see eqn. 3.1) 
one has to find the infimum of 



Since both in Navier-Stokes and MHD the value of ft-min = 1/9 is identical, the 
degree of intermittency is determined by ^'(/imin)- This is both valid for the Eulerian 
as well as for the Lagrangian model which guarantees the monotonicity property. 

Our conclusion from this numerical observation is that the geometry of the most 
singular structures (vortex tubes and current sheets) is not the right quantity to 



hp + 3- D{h) 
l-h 



so that the asymptotic behavior is given by 

+ 3 - D{hjnin) 



P > 1 . 



'inin 
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Table 4. Relative ESS-exponents calculated with respect to the structure function of 

order 2. 




Figure 6. Trajectories with high acceleration in NS (left) and MHD (right, blue denotes 

current density and red vorticity). 

determine the degree of Lagrangian intermittency, but it is more important to look 
at the tracer dynamics. In Figure 6 particle trajectories with high acceleration near 
singular structures are shown. Here the isosurfaces belong to a fixed point in time. 
Important is that in the MHD case the trajectories near the sheet structures are 
smooth. Thus contrary to the Eulerian point of view where the sheets are responsi- 
ble for producing intermittency, they do not contribute significantly to Lagrangian 
intermittency. Large changes with high acceleration occur at the ends of the sheet 
structures. Thus a naive translation of Eulerian to Lagrangian structures is not 
possible. A more detailed investigation of the relation between Lagrangian inter- 
mittency and the small-scale structure of dissipation will be presented elsewhere. 
Using the multifractal approach we also compared the PDFs of velocity and mag- 
netic field increments of the order of the Kolmogorov time with the multifractal 
prediction which is depicted in Figure 7, again assuming the validity of (4.1) for 
the velocity and magnetic field. Here, we have chosen h^i^ = 0.16 to obtain the 
best agreement between the model and prediction. The agreement is not as perfect 
as in the Navier-Stokes case. Here, only the exponential tail could be well described 
by the multifractal model. 



5. Conclusions and open questions 

The presented Navier-Stokes simulations show good agreement with recent experi- 
ments but deviate from predictions of a multifractal model. An observation which 
also could not be described by multifractal modeling is that Lagrangian Navier- 
Stokes intermittency is stronger than in the MHD case whereas the situation is 
reversed for Eulerian statistics. The present situation is depicted in Figure 8. It 
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Figure 7. Comparison of PDFs of velocity(+) and magnetic field(x) increments for small 
r: magnetic field increment are broader. The continuous line corresponds to the multifrac- 
tal approach. 

shows that the multifractal prediction for Lagrangian Navier-Stokes turbulence fits 
well to the simulations of Lagrangian MHD turbulence and vice versa. This again 
shows that it is not easily possible to relate the geometry of the most dissipative 
structures to the strong acceleration events. 
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Figure 8. Comparison of measured scaling exponents to the multifractal prediction 

A second and related critical issue in the Lagrangian treatment is the validity 
of the Kolmogorov like relation I ~ tui which connects Eulerian and Lagrangian 
quantities. At least in the neighborhood of strongly dissipating structures (tubes 
and sheets) this relation has to be altered to I ~ tuq where uq is the mean flow 
produced by the vortex. A similar reasoning was given in [23] to calculate the basic 
mechanism for obtaining exponential tails in the PDF of velocity increments. 

The partial success of the PDF modeling with the multifractal ansatz is mainly 
attributed to the freedom of choosing three free parameters: i) one parameter hidden 
in the normalization of the Reynolds number, ii) the value of /imin and finally iii) 
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the amplitude. If one fixes ftmin to the She-Leveque value hmin = 1/9 then this 
modehng is not able to reproduce the shape of the PDF. 

Thus a deeper understanding of the connection between the geometry of dis- 
sipativc structures and high acceleration events is necessary to correctly model 
Lagrangian intermittency. 
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